1. Field of the Invention
The present invention relates to digital coding of multi-dimensional representations of shapes and objects especially to the representations of three-dimensional shapes or objects and how to acquire such shapes from a plurality of sections through these shapes or objects. The present invention may find advantageous use in the generation, storing and recall of surface representations of multi-dimensional objects, in the transmission of these representations over bandwidth limited transmission channels, in view-dependent decoding or visualization scenarios to adapt the complexity of the scene depending on the processing power in order to satisfy the Quality of Service (QoS), and in editing and animation operations.
2. Description of the Related Technology
Modern medicine cannot be envisaged without the aid of digital imaging. Daily, large amounts of images are obtained from two-dimensional (2D), three-dimensional (3D), four-dimensional (4D—time/space), and multidimensional acquisition devices (i.e. X-ray computed tomography (CT), magnetic resonance imaging (MRI), confocal microscopy, ultrasound imaging, single photon emission computed tomography (SPECT) and positron emission tomography (PET)). While CT and MR acquisition devices still represent a big investment, ultrasound equipment is more accessible for home practitioners and will probably become more and more popular. All these digital images need to be efficiently stored, exchanged and processed. Therefore recently high performance lossless and lossy compression (decompression) algorithms have been developed. Most of these algorithms are based on wavelets and have been tuned in terms of speed and quality with respect to compression ratio.
With these algorithms, compressed images can be exchanged more efficiently, i.e. more data can be available for the same storage capacity, and faster transfer can be achieved via a network channel. Compression schemes, which produce an embedded stream that supports progressive refinement of the decompressed image up to the lossless stage, are well suited for telematics applications that require fast tele-browsing of large image databases while retaining the option of lossless transmission, visualization and archiving of images. This approach permits, in a first step to display only low resolution or low quality images (by sending only a part of the coded stream) in order to get an overview of the contents, and then to request in a second step only the relevant information (i.e. certain slices from a 3D data set or regions of interest) at full image size or full quality by sending the remaining data in the stream without duplication. A tele-medicine environment that enables such a scenario may for example be set up.
The current and future visualization standards go beyond 2D visualization (i.e. orthogonal or arbitrary cross-sections through the 3D data set): they foresee the ability to display 3D models of the objects contained in the images. 3D visualization methods can be divided into two categories: volume rendering approaches and polygonal surface visualization methods. Although the latter methods require an extra surface extraction step, in addition to the segmentation of the objects, they certainly offer a number of benefits:    A better visual quality can be obtained from a shaded surface than from a volume rendering, especially when zooming in for detail.    Shading the surface's triangles is faster than volume rendering when using specialized hardware, making surface rendering more suited for interactive visualization. As an example, graphic boards available for workstations have affordable prices, and currently have the potential to display more than 30 million triangles per second by integrating on the same board engines to compute the geometric transform of the mesh vertices, the lighting calculations and the rasterization.    The meshes of the objects can be exploited in finite element analysis and structure analysis, or inside a CAD package for design, animation and manufacturing purposes. A finite element analysis has been applied for example in M. A. Perzl, H. Schulz, H. G. Paretzke, K. H. Englmeier and J. Heyder, “Reconstruction of the Lung Geometry for the Simulation of Aerosol Transport”, Journal of Aerosol Medicine, 9(3):409-418, 1996, where the mesh representation of a dog's lungs has been used for the simulation of aerosol transport.    A number of mesh editing and post-processing techniques can be applied, such as: simplification, smoothing, intersection or union of meshes, removal of certain parts, etc.    In the medical field, surgery planning is a hot topic. For instance for patients suffering from a brain tumor, usually before surgical intervention several 3D data sets are acquired with different parameters and/or different equipments (MR, CT, . . . ), in order to enhance the geometrical accuracy, visual quality and the automatic segmentation of relevant image structures, such as: tumors, blood vessels, tissue, bones, etc. These data sets have to be co-registered first, so that each position in a certain set can be related to a corresponding position in another set. Since different structures might be visible in different data sets, one has to be able to display them as overlapping images (e.g. displaying the contours of the arteries from a data set on top of the images of another data set) or to generate a 3D view of the brain structures of interest. By first extracting the meshes for the objects of interest from different data sets, and by loading all these meshes together thereafter, one has the ability to visualize in 3D only the relevant information and to plan and simulate the surgical procedure. Also several examples in prosthesis design and placement are currently making use of surface extraction from measured image data sets, mesh generation and visualization (e.g. in dentistry, artificial joints, . . . )    Progressive transmission or refinement of meshes involves coding the mesh description and its vertices in such a way that several resolutions of the same objects can be retrieved from the resulting stream without having redundant information (i.e. each vertex is available only once in the stream and all the vertices from the stream belong to the final resolution object). Such a coding scheme has several applications. For instance, in real-time animations, the complexity of the objects can be adapted depending on the processing power. For some objects not all details are visible. In some situations (e.g. when zooming out image scenes or when dealing with objects located far away from the viewer in a scene with a considerable perspective angle), objects are displayed with a scale factor significantly less than 100%, so that it makes sense to display a coarser (less resource consuming) version of the object's mesh. Progressive transmission of meshes over a network has the benefit that the user at the receiver site does not have to wait until the complete mesh description is downloaded. He/she can start visualizing the object as soon as enough data are available for a coarse level display, and he/she can watch how the model improves while receiving the rest, or he/she may just decide to stop the transfer when he/she sees enough detail.
Since different authors in the field of the present invention use their own notations, the notations used in the present invention (discrete geometry) are given hereinafter for sake of clarity.
First, the notion of a discrete space and a digital image will be defined, and the difference between a digital and a binary image will be shown. Next, iso-contouring is defined and its relation to a binary image is expressed. A digital image consists of voxels, which can be classified either as belonging to the background or to the foreground. The adjacency relation between voxels is introduced and the impact on iso-surface extraction is described.
Let ψD=[1, . . . , δ1]×[1, . . . , δ2]×. . . ×[1, . . . , δn] be an n-dimentional discrete space of size δ1, δ2, . . . , δn. Such a discrete space is called an n-dimensional bounded grid and abbreviated as n-grid. A digital image I (or short image) is an application ψD→D(I) where D(I) is the value domain of the image, obtained from sampling a scalar field σ(x) over the n-grid xε″. Typically D(I)={0,1 for a black and white image (also called binary image), D(I)=[0, . . . , 255] for 8-bit gray level images, and D(I)=[0, . . . , 255]×[0, . . . , 255]×[0, . . . , 255] for 24-bit color images. An n-voxel v is an element of ψD, and the value of v in the image I is defined as I(v).
A common strategy for visualizing a scalar field {overscore (ω)}(x) over an n-grid is to render a (n−1)-dimensional surface satisfying the condition Φ(x)=const for xε″. The visualization technique is called iso-contouring. An iso-contour is defined as C(τ): {x|Φ(x)−τ=0} with τ representing a specified threshold value and Φ(x) representing a computed feature for the n-voxel vx. Such a feature could be the scalar field value σ(x) or any other property extracted from the neighborhood of the n-voxel vx.
For a specified threshold τ, any digital image generally defined as I: D(I)=[a, . . . , b] can be converted to a binary image defined as I′: D(I′)={0,1}, by classifying the n-voxels of the n-grid into two groups: one containing all n-voxels having Φ(x) above or equal to τ, and the second containing the remaining n-voxels having Φ(x) below τ.
Since any digital image can be reduced to a binary image, or to several ones for different computed features Φi(x), the description of the present invention can refer binary images without loosing generality. All the following definitions are valid within a binary picture:    The inverse picture Ī of an image I can be defined as Ī(v)=1−I(v) with vεψD.    The background N(I) of an image can be defined as N(I)={vεψD|I(v)=0}, and the foreground U(I) of an image as: U(I)={vεψD|I(v)=1}.    An object O is a subset of U(I).
Background n-voxels are referred to as o-voxels, respectively foreground n-voxels as •-voxels.
A voxel can be classified as belonging to the object, and located inside the surface, marked as a black dot or •, and referred to as •-voxel, or it can be classified as belonging to the background, and located outside the surface, marked as a white dot or o, and referred to as o-voxel. This is called the state of the voxel.
For the particular case of a 2-dimensional (2D) image, v is called a pixel, and is identified by its two coordinates (i,j)εψD. A pixel can be displayed as a square of unit size. Similar for a 3-dimensional (3D) volume, v is called a voxel and is identified by its three coordinates (i,j, k)εψD.
A voxel can be visualized as a cube of unit size, having 6 faces that are called upper, lower, right, left, front and back, and each face has 4 edges. One is able to slice a 3D volume in orthogonal planes, thus obtaining sequences of 2D images. Since a 2D image consists of pixels, the analogy can be made that each face of a voxel represents a pixel in an orthogonal slice through the volume, as shown in FIG. 1. In FIG. 1(a) an image is shown, containing both background pixels (gray pixels 1) and foreground pixels (white pixels 2). The same image is shown in FIG. 1b as a slice through the volume and containing voxels, gray voxels 3 and white voxels 4. A pixel can be viewed as a voxel's face.
A 3D volume can be created artificially, containing experimental or simulated data for example, or it can be the result of a stack of planar cross-sectional slices produced by an image acquisition device (e.g. MR, CT, confocal microscope, ultrasound, etc.). The pixels within the slices usually have the same resolution in both directions. The slices, however, in many cases are acquired at a different resolution than the pixels within the slice, which results in an anisotropic voxel. For anisotropic volumes, some interpolation within the data set can be done in an extra pre-processing step to achieve the isotropy, or this problem can be solved at display time by applying different scaling factors, which are proportional to the voxel's resolution in {x,y,z}, for the 3D surface vertices.
Two •-voxels are κ-adjacent for κ={6,18,26} if their coordinates differ of ±1:    κ=6: exactly one coordinate (face adjacency) is different, as shown in FIG. 2(a),    κ=18: one or two coordinates (face or edge adjacency) are different, as shown in FIG. 2(b),    κ=26: one, two or three coordinates (face, edge or vertex adjacency) are different, as shown in FIG. 2(c).
Two •-voxels are said to be strictly κ-adjacent if they are adjacent only for this κ. Two κ-adjacent voxels u and v are denoted κ(u,v) and κ is called an adjacency relation. In order to define objects, the transitive closure of this relation, called the connectedness, is used: two •-voxels v1 and vn are κ-connected (in the set ψD) if there exists a sequence of •-voxels v1, v2, . . . , vn such that for any 1≦i<n, vi is κ-adjacent to vi+1. An object O is κ-connected if all pairs of voxels of O are κ-connected.
If κ is specified for a •-voxel, then a similar adjacency relation λ and connectiveness is available for a o-voxel (a description of connectedness couples (κ,λ) is given below).
The κ-neighborhood of a voxel v is a set of voxels that are κ-adjacent with voxel v. Based on this definition, a more general classification of the voxels can be stated. They can be classified into three groups:    Voxels are belonging to the background when their κ-neighborhood contains only o-voxels.    Voxels are belonging to the object when their κ-neighborhood contains only •-voxels.    Voxels are belonging to the border when their κ-neighborhood contains both •-voxels and o-voxels.
The oldest methods to construct polygonal meshes, for the surface representation of a three-dimensional solid, are contour based approaches, which try to connect contours in adjacent 2D slices to form triangular meshes. Usually a lot of user interaction is needed for specifying, in ambiguous cases (i.e. branching), the contours that have to be connected.
These algorithms face three major problems:    The correspondence issue involves finding the correct connections between the contours of adjacent slices.    The tiling problem has two related issues. The first is how to accomplish optimal tiling in terms of certain metrics such as surface area and enclosed volume. The second one is the topological correctness of the tiling, and the detection and tiling of dissimilar portions of contours.    The branching problem occurs when a contour in one slice can correspond to more than one contour in an adjacent slice.
In E. Keppel, “Approximating complex surfaces by triangulation of contour lines”, IBM Journal of Research and Development, 19(1):2-11, January 1975, Keppel provides a solution for building the mesh by triangulation based on randomly distributed points along contour lines. The number of chosen contour points depends on the desired accuracy. The combinatorial problem T of finding the best arrangement of triangles between two contour lines consisting of n and m points can be computed as: T(m,n)=[(m−1)+(n−1)]!/[(m−1)!(n−1)!], which in case of n=m=12, for example, provides 107 combinations and thus different surface shapes.
Keppel's strategy for finding the optimal arrangement of triangles approximating the unknown surface is to decompose the given set of contour points into subsets, according to their position and the shape of the contour, convex or concave, such that a proper objective function may be stated separately for each subset. The objective function is formulated as: The triangulation which maximizes the volumes of the polyhedron gives the optimal approximation of the surface provided by a pair of closed convex contour lines. The objective function can be extended to concave contour point sets as well. The optimal triangulation is found by using classical methods of graph theory. The main unsolved problem is the ambiguity of how to triangulate when branching structures cause splitting of one contour into several ones in the neighboring slice.
In H. Fuchs, Z. M. Kedem, and S. P. Uselton, “Optimal Surface Reconstruction from Planar Contours”, Communications of the ACM, 20(10):693-702, October 1977, Fuchs et al. determine the optimal surface for each pair of consecutive contours by reducing the problem of finding certain minimum cost cycles in a directed toroidal graph. Their algorithm has a time complexity of O(n2 log n), where n is the total number of vertices on the contours bounding the triangles.
In D. Meyers, S. Skinner, and K. Sloan, “Surfaces from contours”, ACM Transactions on Graphics, 11(3):228-258, July 1992, Meyers et al. extend previous work and focus on the correspondence problem (finding the correct connections between the contours of adjacent slices) and on the branching problem. They describe a new approach for solving the correspondence problem using a Minimum Spanning Tree generated from the contours. They approximate the contours by ellipses and then assemble them into cylinders to determine the correspondence. Their branching method handles cases where m contours in one section merge into n contours in an adjacent section. Another improvement is the reduction of the user interventions required to determine where the cross-links between the contours should be. They also discuss a scheme to tile canyons between contours. While their branching method works reasonably well for identifying the openings of canyons between branching contours, it does not handle cases where a contour is enclosed within another and the object interior is in between.
In C. Bajaj, E. Coyle, and K. Lin., “Arbitrary Topology Shape Reconstruction from Planar Cross Sections”, Graphical Models and Image Processing, 58(6):524-543, 1996, Bajaj et al. provide a good survey of previous approaches. In their approach, Bajaj et al. tackle all three problems (correspondence, tiling and branching) simultaneously by imposing three constraints on the reconstructed surface and then deriving precise correspondence and tiling rules from these constraints. In their approach, Bajaj et al. use the 2D marching cubes algorithm as described in W. E. Lorensen and H. E. Cline, “Marching Cubes: A high resolution 3D surface construction algorithm”, in M. C. Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21, pages 163-169, 1987, to generate contour segments from an image slice.
Bajaj et al. give an example of the correspondence problem, shown in FIG. 3. It is to be noticed that four different joint topologies as shown in FIG. 3(b) to (e) may result from the same cross section as in FIG. 3(a). If the distance between slices is large, a priori knowledge or global information is required to determine the correct correspondence.
Tiling means using slice chords to triangulate the strip lying between contours of two adjacent slices into tiling triangles 5, as shown in FIG. 4. A slice chord 6 connects a vertex 7 of a given contour 8 to a vertex 7′ of the contour in an adjacent slice 8′. Each tiling triangle 5 consists of exactly two slice chords 6 and one contour segment.
An example with dissimilar contours 9, 10 located on two adjacent slices is shown in FIG. 5(a). FIG. 5(b) shows a tiling in which all vertices of the top contour 9 tile to vertices of the bottom contour 10. In the vertical cross section shown in FIG. 5(c), which is a cross section passing through the point P of the solid in FIG. 5(b) of a reconstruction as in FIG. 5(b), it is noticed that the scalar data (value) along the vertical line L flips (from outside the surface to inside, or vice versa) twice between two adjacent slices because the surface is intersected twice by L. This is an unlikely topology, especially when the distance between the two slices is small. Another tiling possibility is shown in FIG. 5(d), with the vertical cross section displayed in FIG. 5(e), in which the vertices of the dissimilar portion of the contour tile to the medial axis 11, 11′ of that dissimilar portion, resulting in a highly likely topology.
Bajaj et al. compare their tiling algorithm with previous approaches, taking an example of two dissimilar contours 9, 10 as shown in FIG. 6, which causes many tiling algorithms to fail. A top view of the two contours 9, 10 is illustrated in FIG. 6(a). The thicker contour is the top contour 9, and small circles represent vertices. The result of the minimum surface area-optimizing algorithm described in D. Meyers, “Reconstruction of Surfaces from Planar Contours”, PhD thesis, University of Washington, 1994, is shown as a wire frame in FIG. 6(b) and the same in FIG. 6(c) with hidden lines removed. The arrow points to the abnormality where triangles intersect non-adjacent triangles. It is noticed that the minimum surface-optimizing algorithm generates non-polyhedral surfaces. The result of the shortest slice cord heuristic algorithm as described in H. N. Christiansen and T. W. Sederberg, “Conversion of complex contour line definitions into polygonal element mosaics”, Computer Grapics, 12:187-192, August 1978, is shown in FIG. 6(d), which fails as well. Even when the tiling result is a polyhedron, it might be physically unlikely, as shown in FIG. 6(e), where there exists a non-self intersecting surface. One polyhedron correct solution obtained by the algorithm of Bajaj et al. is shown in FIG. 6(f).
A branching problem is shown in FIG. 7, where one contour C3 of slice S2 branches into two contours C1 and C2 of slice S1 (see FIG. 7(a)). FIGS. 7(b) and 7(c) show that a curve L or a point is added between the two slices S1, S2 to model the valley or saddle point formed by the branching. This method is used by Bajaj et al. and is described in C. Bajaj, E. Coyle, and K. Lin., “Arbitrary Topology Shape Reconstruction from Planar Cross Sections”, Graphical Models and Image Processing, 58(6):524-543, 1996. In H. N. Christiansen and T. W. Sederberg, “Conversion of complex contour line definitions into polygonal element mosaics”, Computer Grapics, 12:187-192, August 1978, Christiansen et al. use the method as shown in FIG. 7(d). This approach works well only in simple branching cases. In D. Meyers, S. Skinner, and K. Sloan, “Surfaces from contours”, ACM Transactions on Graphics, 11 (3):228-258, July 1992, Meyers et al. use the scheme as shown in FIG. 7(e), and afterwards they improve the horizontal triangle problem associated with this approach by feeding the triangulation mesh into a surface-fitting program to regenerate the surface.
The handling of branching according to Bajaj et al. is shown in FIG. 8. Branching regions that are not complex, as in the top view shown in FIG. 8(a), can be tiled by adding a vertex to model the saddle point, as shown in FIG. 8(b). The thicker line segments in FIG. 8(a) are top contours 9, 9′, and small circles represent vertices. Alternatively, for example for the branching contours 9, 9′, 9″, 10 shown in top view in FIG. 8(c) (again the thicker line segments are top contours 9, 9′, 9″, and small circles represent vertices), one could use line segments 13 to model canyons between the contours, as shown in FIG. 8(d).
The constraints imposed by Bajaj et al., ensure that the regions tiled by these rules have a natural appearance and the surfaces produced correspond well with expected physical models. For regions that cannot be tiled by these rules without violating one or more constraints, the medial axis of the region is considered and the tiling is done correspondingly.
The surface reconstruction criteria defined by Bajaj et al. are the following:    Criterion 1: The reconstructed surface and solid regions form piecewise closed surfaces of polyhedra.    Criterion 2: Any vertical line (a line perpendicular to a slice) between two slices intersects the reconstructed surface at zero points, at one point, or along one line segment.    Criterion 3: Resampling of the reconstructed surface on the slice should produce the original contours.
These criteria are used to derive correspondence and tiling rules. The correspondence rules are local (data in adjacent slices are used to determine the correspondence between contours). The tiling rules prohibit those undesired tilings or nonsensical surfaces, and allow detection of branching regions and dissimilar portions of contours. In conjunction with the rules, the method combines a multipass tiling algorithm that first constructs tilings for all regions not violating any of the tiling rules, and in a second and third pass, processes regions that violate these rules (holes, branching regions and dissimilar portions of contours), by tiling to their medial axes. Topologies that violate criterion 2 due to under-sampling of the original data might generate distorted results, as shown in FIG. 9, where examples of topologies are shown that cannot be processed because they violate criterion 2 by having two points of intersection with a line L perpendicular to the slice. An example of surface reconstruction of the brain hemisphere produced by the algorithm of Bajaj et al. is shown in FIG. 10. It is generated from a set of contour data that has been manually traced from 52 MRI image slices. The approach produces significantly fewer triangles than the marching cubes approach as described in W. E. Lorensen and H. E. Cline, “Marching Cubes: A high resolution 3D surface construction algorithm”, in M. C. Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21, pages 163-169, 1987, and as described in U.S. Pat. No. 4,710,876, but is much slower. FIG. 10(a) represents a Gouraud shading and FIG. 10(b) represents a wire frame of a reconstructed brain hemisphere.
The current standard technique for extracting iso-surfaces (Σ) from a 3D grid of data is the Marching-Cubes algorithm (MC), revealed in the patents U.S. Pat. No. 4,710,876 and U.S. Pat. No. 4,729,098, and the related paper W. E. Lorensen and H. E. Cline, “Marching Cubes: A high resolution 3D surface construction algorithm”, in M. C. Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21, pages 163-169, 1987.
MC is based on the idea of combining interpolation, segmentation, and surface construction in one single step, but this implicit segmentation is sometimes very difficult, if not impossible, e.g. with noisy data.
At the time of the MC algorithm introduction, existing Σ extraction algorithms lacked in detail and sometimes introduced artifacts, as discussed in E. Keppel, “Approximating complex surfaces by triangulation of contour lines”, IBM Journal of Research and Development, 19(1):2-11, January 1975, and in H. Fuchs, Z. M. Kedem, and S. P. Uselton, “Optimal Surface Reconstruction from Planar Contours”, Communications of the ACM, 20(10):693-702, October 1977. MC creates a polygonal representation of constant density surfaces from a 3D grid data, producing models with “unprecedented detail”, as described by the authors. The original MC algorithm has been modified by various fixes for its failure to construct a true separating surface. This has been described in H. H. Baker, “Building surfaces of evolution: The weaving wall”, Int. J. Comp. Vision 3 (1989), 51-71; in A. Van Gelder and J. Wilhelms, “Topological Considerations in Iso-surface Generation”, to appear in ACM Transactions on Graphics; in A. D. Kalvin, “Segmentation and surface-based modeling of objects in three-dimensional biomedical images”, PhD thesis, New York University, 1991; in G. M. Nielson and B. Hermann, “The asymptotic decider: resolving the ambiguity in marching cubes”, Proceedings of Visualization '91 (1991), 83-90; in P. Ning and J. Bloomenthal, “An Evaluation of Implicit Surface Tilers”, IEEE Computer Graphics & Applications, 13(6):33-41, November 1993; and in B. Wyvill and D. Jevans, “Table Driven Polygonization”, SIG-GRAPH 1990 course notes 23, 7-1-7-6.
The “original” MC algorithm as described in U.S. Pat. No. 4,710,876 and U.S. Pat. No. 4,729,098 is based on the tabulation of 256 different configurations, thus it falls into the category of methods that are based on a lookup table to decide how to connect the iso-surface vertices. It allows a local iso-surface extraction inside an 8-cube defined by eight voxels. The global iso-surface is the union of these small pieces of surface computed on each set of eight connected voxels of the image.
An 8-cube as shown in FIG. 11(a), is a cube 14 lying between two image slices 15, 16 and consisting of 8 neighboring voxels (4 in each slice). A new grid is formed by identifying the image voxels as discrete 3D points. For a 3D image I with size M×N×P, there are (M−1)(N−1)(P−1) 8-cubes in I. The voxels, or corners, of the 8-cube 14 are either located inside (•-voxel) or outside (o-voxel) of the object 17. Whenever there are two corners with different states (inside and outside) along an edge of the 8-cube 14, there must be an intersection with the surface of the object 17. Since there are 8 corners (voxels) in each 8-cube 14, and two states, inside and outside, there are 28=256 ways a surface can intersect the 8-cube 14. Triangulating these 256 cases is possible but tedious and error-prone. That is why these cases are further reduced to 15, first by eliminating complementary cases (interchanging •-voxels with o-voxels and vice versa) and second by omitting symmetrical cases (obtained by rotating the 8-cube). FIG. 12 shows the basic 8-cube configurations for MC.
In order to speed-up the identification of triangles that are built inside an 8-cube, a lookup table with 256 entries has been generated from these 15 cases. The lookup table contains a list of 8-cube edge intersections with a description stating how these intersections should be connected in order to build the triangles. For each case, an index into the lookup is computed using the corner numbering of the 8-cube as shown in FIG. 11(b), by associating each corner with one bit in an 8-digit number, and assigning a one to each •-voxel and a zero to each o-voxel belonging to the 8-cube. For the example shown in FIG. 11(b), which corresponds to case 5 of FIG. 12, giving rise to three triangles, the computed index=01100100 binary (100 decimal), and the connectivity of the triangles, described in terms of edges intersected by Σ, is {e3e5e7,e2e3e5,e2e5e10}.
The position of each triangle vertex, where the iso-surfaces Σ intersects the 8-cube's edge, is found by linear interpolation of the density (e.g. gray value) between the edge corners. The computed offset from the •-voxel towards the o-voxel for a specified threshold τ is equal to:                     offset        =                              (                                          D                ⁡                                  (                                      •                    -                    voxel                                    )                                            -              τ                        )                                (                                          D                ⁡                                  (                                      •                    -                    voxel                                    )                                            -                              D                ⁡                                  (                                      •                    -                    voxel                                    )                                                      )                                              (eq.  1)            where D is the density value at the corner voxel. The offset is in the range=[0 . . . 1].
For visualization purposes, a unit normal is computed for each triangle vertex using the gradient information available in the 3D data. First, the gradient vector at the corners of the 8-cube is estimated and then the vectors are linearly interpolated at the point of intersection of the triangle vertices with the 8-cube edges. The gradient at the 8-cube corner (x,y,z), is estimated using central differences along the three coordinate axes by:Gx(x,y,z)=[D(x+1,y,z)−D(x−1,y,z)]/2ΔxGy(x,y,z)=[D(x,y+1,z)−D(x,y−1,z)]/2ΔyGz(x,y,z)=[D(x,y,z+1)−D(x,y,z−1)]/2Δzwhere D(x,y,z) is the density at voxel (x,y,z) and {Δx,Δy,Δz} are the lengths of the 8-cube edges.
The MC algorithm has been intensively used for rendering purposes but has shown its limitations in other applications, because the extracted iso-surface is in general not a “simple” surface. Some authors have noticed that the iso-surface does not always have the same topology as an underlying continuous surface defining the image. Many authors have contributed to solve this problem, but they have often provided an empirical solution or a visual justification.
There are some ambiguous cases that might lead to the generation of holes, i.e. inconsistencies in the produced mesh. In case the resulting mesh is aimed at direct rendering, inconsistencies might not be an issue, as long as one is not zooming in for details. In a different context, for instance if the goal is to use the mesh as a CAD model or for a finite element analysis, as in M. A. Perzl, H. Schultz, H. G. Paretzke, K. H. Engimeier, and J. Heyder, “Reconstruction of the Lung Geometry for the Simulation of Aerosol Transport”, Journal of Aerosol Medicine, 9(3):409-418, 1996, one might encounter serious problems.
The ambiguous cases occur on any 8-cube face that has adjacent vertices with different states, while diagonal vertices are in the same state. There are six of these cases as shown in FIG. 12: cases 3, 6, 7, 9, 10 and 13. An example of inconsistent reconstruction of the local surface and of a correct solution are shown in FIG. 13. In the left image, an inconsistent choice is shown. Case 6 of FIG. 12 is applied for the complementary configuration that, combined with the triangles generated by case 3, will result in holes, as triangles generated by case 3 do not match the ones generated by case 6 (actually the complement of case 6). A correct triangulation is shown in the right image where case 6 has been modified, a configuration that is not available in the 15 basic cases.
A first fix to the MC algorithm, is constructing iso-surfaces from CT data, as described in A. Wallin, “Constructing isosurfaces from CT data”, IEEE Computer Graphics and Applications, 11(6)28-33, November 1991.
Similar to MC, this method uses the 8-cube approach to extract a local surface. The main difference is that the configuration table is built for a 4-face instead of an 8-cube. The ambiguous cases found in MC are solved by using the density value (interpolated) at the center of the 4-face to determine if that point lies inside or outside the Σ, as described in G. Wyvill, C. McPheeters, and B. Wyvill, “Data structures for soft objects”, The Visual Computer, 2:227-234, 1986.
Some new concepts are introduced. A polygon is a connected set of edges that forms an iso-surface. The ordering of a polygon is the path taken by the vertices describing it (v0→v1→v2→ . . . →v0). Two neighboring polygons are coherently ordered if the path of an edge shared by the two is opposite. A surface is coherently connected if all edges occur twice in opposite directions (two complementary edges) and if no polygon touches another polygon except at the common edge.
The algorithm guarantees that the obtained polygons are coherently ordered and connected, with no polygon occurring more than once, and that each surface is complete, that is, no holes due to surface generation errors occur. The generated Σ consists of either the minimum number of (non-planar) polygons with 3 to 12 vertices or a minimum number of triangles.
The surface generation process has two phases that are: the edge generation and connection, and polygon generation.
In the first step, for each 8-cube three 4-faces (for example the front, right and top 4-faces) are considered and the edges intersected by Σ are determined. The other three of the six 4-faces (for example the back, left and bottom 4-faces) will be examined when the adjacent 8-cube are considered. With four voxels at the corners of the 4-face, 16 cases of edge configurations must be considered, as shown in FIG. 14. White dots represent the o-voxels, black dots represent the •-voxels, and gray dots represent the vertices of the Σ. For configurations 5 and 10, either the dashed or the solid lines connecting grey dots indicate which edges are used, depending on the (interpolated) value at the center of the face. In contrast to MC, the cases are not reduced by symmetry or complementariness.
When one intersection occurs, the algorithm produces two edges (case 1 to 4, 6 to 9, and 11 to 14), one directed in each direction. Both of these edges must be kept in memory to produce a coherent surface later. When two intersections occur (as in case 5 and 10), the algorithm can use the (interpolated) value at the center of the face to determine if that point lies inside or outside the iso-surface, as described in G. Wyvill, C. McPheeters, and B. Wyvill, “Data structures for soft objects”, The Visual Computer, 2:227-234, 1986. The generated surface differs slightly depending on whether the face's center is considered or not. Either way the obtained surface is valid and coherent.
Similar to MC, the position of the triangle vertex, where the iso-surfaces Σ intersect the 4-face's edge, is linearly interpolated between the edge corners.
In the second step, the scanning of the edges list to obtain a correct surface is more complex than the generation of edges. The possible polygon configurations that can be detected are triangles (one case), rectangles (two cases), pentagons (one case), hexagons (three cases), heptagons (one case), octagons (two cases), nonagons (one case), and 12 edges (two cases). In general, the polygons formed by connecting the edges are non-planar and the triangulation problem is not trivial since the number of possible (although not all valid) configurations is:   N  =                    2                  n          -          2                                      (                      n            -            1                    )                !              ⁢                  ∏                  k          =          1                          n          -          3                    ⁢              (                              2            ⁢            k                    +          1                )            where n is the number of vertices in the polygon. For n=4 the number is 2, for n=6, it is 14, but for n=12 it is 16796.
A second fix to the MC algorithm is skeleton climbing (SC), as described in T. Poston, H. T. Nguyen, P.-A. Heng, and T.-T. Wong, “Skeleton Climbing: Fast Isosurfaces with Fewer Triangles”, Proceedings of Pacific Graphics'97, pages 117-126, Seoul, Korea, October, 1997, and in T. Poston, T.-T. Wong and P.-A. Heng, “Multiresolution Isosurface Extraction with Adaptive Skeleton Climbing”, Computer Graphics Forum, 17(3):137-148, September 1998.
Similar to the Marching Cubes algorithm (MC), SC uses the 8-cube approach to visit the data volume and to build triangulated iso-surfaces (Σ) in 3D grid data. It is based on a 4-step approach to construct the surface:    First, surface vertices (0-skeleton of the Σ) are chosen just as MC does, i.e. one on each cube edge joining a o to a •, positioned by linear interpolation. An edge is called occupied if it contains a surface vertex. As an auxiliary for the triangle building step, thick edges 18 as shown in FIG. 15 are introduced and defined in the following way: in z-direction the ones with x and y even, in y-direction the ones with x odd and z even, and in x-direction the ones with y and z odd.    Rather than to look directly at the configuration of the cubes like MC does, triangle sides (1-skeleton of the Σ) are constructed in the second step between vertices that share a cube's face. A simple rule, “do not cross o—o diagonals”, which yields five different cases for the face as shown in FIG. 15, is specified. Gray dots represent the vertices located between a o- and a •-voxel, which are connected by oriented edges. One can notice that SC does not reduce the number of cases by symmetry and complementariness, so that in configurations similar to case 4 in FIG. 15 mismatching edges can be avoided, a problem that the original MC algorithm faces, resulting in holes in the surface. Each side has an orientation such that the neighboring •'s stay on the left side.    In a third step, polygons are retrieved in each cube of the cubical grid and split into triangles (2-skeleton of the Σ) according to the following rule: whenever the edge E of a cube is occupied, the set of sides is shrunk by merging the two sides that end on E, and a triangle is generated with the orientation given by these two sides, as shown in FIG. 16.    After the third step, the triangulated surface is fully defined. The fourth step is only supposed to reduce the number of triangles: for each occupied thick edge 18, the four triangles that meet on it, are replaced by two triangles with the same boundary, as shown in FIG. 17. The four triangles in FIG. 17(a), resp. FIG. 17(c), that meet on the thick edge 18 are replaced by two triangles with the same boundary in FIG. 17(b), resp. FIG. 17(d). FIGS. 17(c) and (d) show the same configuration and procedure as in the FIGS. 17(a) resp. (b), but in the context of the full mesh. The method allows still to obtain truly separated surfaces without gaps while reducing the number of triangles by 25%. The time overhead introduced by this step is 4% on average, yielding globally a time comparable with MC, but that is however significantly smaller than the time needed by current mesh decimation algorithms.
The authors emphasize that they do not reduce the cases by symmetry and complementarity in the SC algorithm as MC does. It is possible to derive 256 configurations from the five face cases shown in FIG. 15, but to make a comparison possible with MC, it suffices to only look at the basic configurations (each one representing equivalent cases by reflection and rotation) in FIG. 18, which shows typical SC step 3 triangle sets and side patterns in cube faces for the “do not cross o—o diagonals” rule. The numbering and views are kept similar to those in W. E. Lorensen and H. E. Cline, “Marching Cubes: A high-resolution 3D surface construction algorithm”, in M. C. Stone, editor, Computer Graphics (SIGGRAPH '87 Proceedings), volume 21, pages 163-169, 1987 (Case 14 in the MC literature, identical to case 11 by reflection in a diagonal plane, is omitted here). It is to be noticed that for the cases 3 versus 3a, respectively 6 versus 6a, and 7 versus 7a, complementation gives a different triangulation pattern, an aspect that MC is not taking into account. The other cases are similar to the ones in MC.
Some results that compare the SC algorithm with MC are shown in FIG. 19 and FIG. 20. FIG. 19 is an example of a knotted torus sampled at 64×64 resolution, whereby FIG. 19(a) is extracted by MC (13968 triangles), FIG. 19(b) is extracted by SC without step 4 (13968 triangles), and FIG. 19(c) is extracted by SC with step 4 (10464 triangles). Without step 4, SC obtains similar results as MC in terms of number of triangles, speed and visual quality, except for the cases where MC fails to correctly approximate Σ. SC with step 4 reduces the number of triangles by 25%, and gives as a result a smoother surface, as also shown in FIG. 20, except for fine details, which are approximated less accurately; see FIG. 19(c) where step 4 preserves the topology, but for an analytic function, some smoothness is lost in the finest part of the tube. FIG. 20 shows surfaces extracted from a 128×128×57 CT scan, with a threshold distinguishing bone from non-bone. Again FIG. 20(a) is extracted by MC (100830 triangles), FIG. 19(b) is extracted by SC without step 4 (100066 triangles), and FIG. 19(c) is extracted by SC with step 4 (75318 triangles). SC with step 4 gives a smoother forehead region, and equally good detail.
As described in T. Poston, T.-T. Wong and P.-A. Heng, “Multiresolution Isosurface Extraction with Adaptive Skeleton Climbing”, Computer Graphics Forum, 17(3):137-148, September 1998, Poston et al. have implemented a new SC algorithm called adaptive skeleton climbing (ASC) that is able to generate directly multi-resolution iso-surfaces from volume data, with between 4 to 25 times fewer triangles than MC, and in comparable processing times. The basic idea is to group voxels first in 1D (segments), then in 2D (rectangles) and finally in 3D (boxes) that are adapted to the geometry of the true iso-surface. In ASC, the triangular mesh is generated from these boxes instead of the cubes as in SC. The coarseness of the generated meshes is controlled by a specified maximum size (N) of the boxes. When N=1, the largest box contains 2×2×2 voxels, like in MC. When a larger block size is used, larger boxes will be generated when the iso-surface can be respected, resulting in larger-triangles.
The proposed on-the-fly triangle reduction approach can generate meshes that are more accurate because it directly makes use of the voxel values in the volume. The algorithm has the advantage that the construction of coarser iso-surfaces requires less computation time, which is contrary to mesh optimization approaches, generally requiring a longer time to generate coarser meshes. Hence, ASC allows the user to obtain rapidly a low-resolution mesh for interactive preview purposes and to delay the decision for the more time-consuming generation of the detailed high-resolution mesh until optimal parameter settings (threshold, viewing angle, . . . ) have been found.
Although the triangles in the meshes may not be optimally reduced, ASC is much faster than post-processing triangle reduction algorithms. Therefore, the coarse meshes ASC produces can be used as the starting point for mesh optimization algorithms, if mesh optimality is the main concern. Similar to SC, ASC does not suffer from the gap-filling problem.
A third fix to the MC algorithm is topologically defined iso-surfaces, as described in J. O. Lachaud, “Extraction de surfaces à partir d'images tridimensionnelles: approches discretes et approche par modèle deformable”, Ph.D. Thesis, July 1998.
In this approach, the authors have focused on modifying the table of configurations for MC, while the algorithm itself remained unchanged. They exploit the underlying discrete topology of voxels, especially their adjacency and connectedness to extract a topologically correct iso-surface from a volumetric image.
Their main contribution is to provide a formal proof of the validity of the generated iso-surface according to the chosen connectedness. Furthermore, they demonstrate the coherence of the iso-surface (closedness, orientability, no singularity, no self-intersection) along with fundamental properties:    For a specified connectedness couple (κ,λ), the obtained iso-surface is opposite (in terms of orientation) of the iso-surface extracted for the opposite image with inverse connectedness couple (λ,κ) (next paragraphs explain the meaning of (κ,λ) and which are the valid connectedness couples).    The •-voxel with a given connectedness is inside the surface and the o-voxel (with another connectedness) is outside the surface.
As shown in FIG. 2, there are three types of adjacencies between image voxels. The adjacency relation can be specified for both •-voxel and o-voxel, call it κ respectively λ.
The orientation of the polygon's edge is important in defining the orientation of the surface. The concept of the orientation of the edges is demonstrated in FIG. 21 for the 16 cases of a 4-face, while the generalization for the 8-cube is shown in FIG. 23. The inner part of the surface (see FIG. 21), displayed in gray, is located on the left side of the oriented edge. Two cases (case 5 and 10) are ambiguous, and the choice of the right connectedness couple (κ,λ) is essential.
A correct solution cannot be obtained for every couple (κ,λ) as shown in FIG. 22. The case of FIG. 22(b) where (κ,λ)=(6,18) or (6,26) and the case of FIG. 22(c) where (κ,λ)=(18,6) or (26,6) generate a coherent iso-surface, while for the case of FIG. 22(a) where (κ,λ)=(6,6), and for the case of FIG. 22(d) where (κ,λ)=(18,18), (18,26), (26,18), or (26,26), an arbitrary choice would provide a locally different surface in the inverse picture. The connectedness couples (6,18), (6,26), (18,6) and (26,6) are called valid couples. It is to be noticed that no distinction can be made between the 18-connectedness and 26-connectedness for a 4-face, but the difference is visible for the 8-cube. The conclusion is that one could build four different tables of configurations for each of the valid couples.
The modified MC configurations for the 14 cases and the valid connectedness couples are shown in FIG. 23. Similar to original MC, cases that can be obtained via rotation or symmetry are not displayed. FIG. 23(a) shows all possible cases obtained for the couples {(κ=6,λ=18), (κ=6 ,λ=26)}. In FIGS. 23(b) and (c) only the cases that differ from FIG. 23(a) for the couple {(κ=18,λ=6) are displayed, with in FIG. 23(c) only cases different from FIG. 23(a) in the way the generated polygon is split into triangles are given. Finally in FIG. 23(d) are shown only the cases different from FIG. 23(b) and FIG. 23(c) for the couple {(κ=26,λ=6).
The result of the modified MC algorithm for different connectedness couples (κ,λ) on a test image is shown in FIG. 24. FIG. 24(a) displays a 3D surface rendering of the individual voxels, whereby each voxel is shown as a box. FIG. 24(b) shows the result obtained for a connectedness couple (κ=26,λ=6), while FIG. 24(c) is obtained for (κ=18,λ=6), and FIG. 24(d) for (κ=6,λ=18) or (κ=6,λ=26).
It is an object of the present invention to provide a method and a system to acquire surface descriptions from a plurality of sections through multi-dimensional representations of objects, especially through three-dimensional representations of objects.
It is a further object of the present invention to provide a compact, optionally multi-scalable, optionally view-dependent, optionally animation-friendly multi-dimensional surface representation method and system.
It is still a further object of the present invention to provide a digital coding of a surface representation which may be transmitted over bandwidth limited communication channels.